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Abstract —The short-term forecasting of real-time locational 
marginal price (LMP) and network congestion is considered from 
a system operator perspective. A new prohahilistlc forecasting 
technique is proposed based on a multiparametric programming 
formulation that partitions the uncertainty parameter space into 
critical regions from which the conditional prohahllity distrihu- 
tlon of the real-time LMP/congestion is obtained. The proposed 
method incorporates load/generation forecast, time varying oper¬ 
ation constraints, and contingency models. By shifting the com¬ 
putation cost associated with multiparametric programs offline, 
the online computation cost is significantly reduced. An online 
simulation technique by generating critical regions dynamically 
is also proposed, which results in several orders of magnitude 
improvement in the computational cost over standard Monte 
Carlo methods. 

Index Terms —Congestion forecast, electricity price forecast, 
locational marginal price (LMP), multiparametric programming, 
probabilistic forecast. 


I. Introduction 

As more renewable resources are integrated into the power 
system, and the transmission system operates closer to its 
capacity, congestion conditions become less predictable and 
locational marginal prices (LMP) more volatile. The increased 
congestion and LMP uncertainties pose significant challenges 
to the operator and market participants, which motivates us 
to consider the problem of short-term forecasting real-time 
LMP/congestion in the presence of generation, demand, and 
operation uncertainties. 

The benefit of accurate LMP and congestion forecasts is 
twofold. For market participants, accurate forecast of real¬ 
time prices is valuable in risk management, bidding strategy 
development, and demand side participation. The forecast 
prices allow market participants to make adjustments in ad¬ 
vance to ensure competitive transactions. For system operators, 
on the other hand, forecast of transmission congestion is 
important for congestion management and system planning. 
European transmission system operators, for instance, use 
Intraday Congestion Forecast (IDCF) to improve real-time 
security assessment a. Intuitively, an LMP forecast should 
elicit generation participation at times of potential shortage 
thus alleviating future congestions. Similarly, an LMP forecast 
can be used for demand response that results in shifting part 
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of the load from peak to valley. An example of applying LMP 
forecast for inter-regional interchange is presented in la. 

Currently, some system operators are providing real-time 
price forecasts. The Electric Reliability Council of Texas 
(ERCOT) im offers a 1-hour ahead real-time LMP forecast, 
updated every 5 minutes. The Alberta Electric System Oper¬ 
ator (AESO) 0 provides two short-term price forecasts with 
prediction horizons of 2 hours and 6 hours, respectively. 

Most LMP forecasting schemes provide only point forecast, 
which gives a single quantity as the prediction. Eor systems 
with high levels of uncertainty, a point forecast is rarely accu¬ 
rate, and impacts of prediction error are difficult to quantify. 
A more attractive alternative is a probabilistic forecast that 
provides a full characterization of the LMP distribution. 

Significant technical challenges exist for probabilistic fore¬ 
casting real-time LMP and congestion. Eirst, reasonably accu¬ 
rate models for real-time dispatch and LMP are needed. Sec¬ 
ond, network operating conditions and uncertainties need to 
be incorporated in real time. Einally, the forecasting algorithm 
needs to be simple and scalable for sufficiently large systems. 

These challenges are daunting for external market partici¬ 
pants who do not have access to network operating conditions 
and confidential information on bids and offers that influence 
LMPs. On the other hand, if it is the system operator providing 
the forecast, as in the case of ERCOT or AESO, the barrier 
to efficient and accurate forecast is lowered. 

A. Summary of Contributions 

In this paper, we consider the real-time LMP and congestion 
forecasting problem from an operator perspective. We focus 
on probabilistic forecasting that, at time t, provides the con¬ 
ditional probability distribution at time t + T of the LMP 
vector and associated congestion, given the system state at 
time t. Here T is referred to as the prediction horizon which 
is considered in the range of T from 1 to 6 hours for short-term 
forecasts. 

The key idea behind the proposed approach is the use of 
multiparametric programming that partitions the uncertainty 
space into critical regions with each region attached to a 
unique LMP and congestion pattern. Thus, the problem of 
probabilistic forecasting reduces to computing the probabilities 
of random parameters falling in the set of critical regions. 
When loads or generations (treated as negative loads) are ran¬ 
dom, their forecasts are incorporated to generate probabilistic 
LMP and congestion forecasts. 

The proposed technique also provides several new fea¬ 
tures not present in existing methods. Eor example, it can 
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incorporate system contingency models and allow system 
constraints to vary with time. The latter feature is relevant 
because network topology and thermal limits may be changed 
in real time by the system operator depending on operating 
conditions. In terms of the generation cost, the proposal can 
be applied to a linear (or piece-wise affine) function and 
a quadratic function. Computationally, the proposed method 
shifts a majority of the computation offline, which significantly 
reduces online computation cost. 

An alternative algorithm that dynamically generates critical 
regions is also proposed. Because load and stochastic gen¬ 
eration processes are physical processes, they are bounded 
and tend to concentrate around the mean trajectory. Their 
realizations thus only fall in a few critical regions instead of 
all over the entire parameter space. By generating the critical 
regions that contain such realizations, the computation cost 
is reduced by several orders of magnitude comparing with 
standard Monte Carlo techniques. 

B. Related Work 

Much of the existing work deals with point forecasts of 
LMP by market participants who do not have access to real¬ 
time operating conditions and confidential offers and bids. For 
these techniques, historical data on LMP, load, and congestions 
drive the forecasting engine. Literature on these techniques 
abounds. See m and references therein. For probabilistic 
forecasting techniques by market participants, see approaches 
in Global Energy Forecasting Competition 2014 Q. 

There are several prior studies on LMP/congestion fore¬ 
casting from the system operator perspective. The proposed 
technique in 0 employs an online Monte Carlo sampling 
technique that, for each Monte Carlo sample path, solves an 
optimal power flow (OPF) problem, which is computationally 
expensive. Monte Carlo technique was also used in ||9| where 
a reduction of the random variable dimension is made using 
a nonhomogeneous Markov chain model based on a partition 
of the system state space. 

A particularly relevant prior work is nni where the authors 
consider the problem of LMP/congestion forecasting from 
the vantage point of an external observer who has access to 
publicly available historical data only. Our work, in contrast, 
considers the forecasting problem from the vantage point of 
a system operator who has access to the system operating 
condition at the time of forecasting. In terms of forecasting 
methodology, the main difference between our approach and 
that in nni lies in the different uses of conditioning in evaluat¬ 
ing the conditional probability distribution of LMP/congestion. 

The authors of cni introduce and exploit the decomposition 
of a multi-dimensional load space into critical regions (called 
system pattern regions) that are estimated using historical 
dat^. The work in ifTOll aims to address the following issue; 
Given a possible future point L in a multi-dimensional load 
space, what is the probability distribution of the estimated 
critical regions that contain L? Since each critical region 
corresponds to a specific LMP/congestion, the technique in 
ifTOll gives a heuristic estimate of the probability distribution 

*The estimated critical regions are therefore random quantities. 


of LMP/congestion by conditioning on load L at some future 
point in time. 

In contrast to cni, our objective is to forecast directly the 
probability distribution of future LMP/congestion, conditional 
on the current system operating point. Because a system 
operator has access to all private and public information about 
system conditions, the critical regions are computed exactly 
via a multiparametric program. This allows us to incorporate 
load and generation forecasts and obtain the (conditional) 
probability distribution of future LMP/congestion directly. 

Several techniques have been proposed to approximate the 
LMP distribution at a future time. In im, a probabilistic 
LMP forecasting approach was proposed based on attaching 
a Gaussian distribution to a point estimate. The advantage is 
that the technique can incorporate various point forecasting 
methods. The disadvantage, on the other hand, is that network 
effects are not easy to incorporate. The authors of ifTSll and 
ifTSll approximate the probabilistic distribution of LMP using 
higher order moments and cumulants. These methods are 
based on representing the probability distribution as an infinite 
series involving moments or cumulants. In practice, computing 
or estimating higher order moments and cumulants are very 
difficult; lower order approximations are necessary. 

C. Organization and Notations 

This paper is organized as follows. Section HI] introduces 
a model of the real-time economic dispatch and the ex-ante 
LMP computation. The key (known) results in multiparametric 
programming that form the basis of the proposed probabilistic 
forecasting approach are discussed in Section |III| Details of 
the proposed techniques are given in Section |IV| and Section 
rvl Numerical results are presented in Section |yT| and followed 
by some concluding remarks in Section IVllI 

The notations used in this paper are standard. For a random 
variable x, its expected value is denoted by x. We use 0 to 
denote an estimate of parameter 0. For a random process yt, 
the forecast of yt+T using all the information available at time 
t is denoted by yt+T\t^ where T is the prediction horizon. We 
use yt+T\t to denote the conditional mean of yt+r given yt 
(and possibly additional information at time t). The notation 
z ~ ^{p,, E) means that z is a Gaussian random vector with 
mean p and covariance matrix S. The indicator function Iyi(a;) 
is one if x € A and zero otherwise. 

Vectors and matrices are in the real field. /„ denotes an 
n X n identity matrix and 1 an n dimensional column vector 
with all ones, where n is the number of buses in the system. 0 
denotes a matrix/vector of zeros with different but compatible 
dimensions. We use superscript “T” to denote the transpose of 
a matrix and the complement of a set J is denoted by J'^. 

11. Real-Time Economic Dispatch and Ex-Ante LMP 

We consider an ex-ante LMP model that arises from a 
real-time economic dispatch. Our model is adopted from the 
stylized model in na that captures the process of com¬ 
puting LMP by independent system operators. Specifically, 
the system operator solves a DC-OPE problem for optimal 
economic generation adjustment that meets load and stochastic 
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generation forecasts for the next dispatch interval and satisfies 
generation and transmission constraints. 

We describe here a simplified real-time economic dispatch 
formulation for the ex-ante LMF0 model. For simplicity, we 
assume that each bus has a traditional generator, a stochastic 
generating unit and a load. The DC-OPF problem at time t is 
given by 

min c^g 
a 

subject to 

V{g + wt+i\t - dt+i\t) = 0 i^t+i) 

L~ < S{g + wt+i\t - dt+i\t) < L+ 

G- <g<G+ 

( 1 ) 

where 

c vector of real-time generation offers; 

g vector of ex-ante dispatch at time f -f 1 ; 

dt+i\t vector of 1 -step ahead load forecast at time f; 

vector of 1 -step ahead forecast of stochastic gen¬ 
eration at time t\ 

S shift factor matrix; 

G^IG~ vector of max/min generation capacities; 

!L~ vector of max/min transmission limits; 

At+i shadow price for the energy balance constraint at 
time f -f 1 ; 

y^t+i shadow prices associated with max/min transmis¬ 
sion constraints at time t -\-l. 

The stochastic generation referenced above can be of any 
form of renewable integrations including renewable energy 
generation, distributed generation, and demand response. Here 
we assume that such stochastic generation is non-dispatchable 
with possible curtailment. The proposed forecasting method 
applies to other cost functions such as a piece-wise affine 
function or a quadratic function. This will be further addressed 
in the following sections 

The LMP value at each bus is the sum of the marginal price 
of generation at the reference bus and the marginal congestion 
price at the location associated with the active transmission 
constraints. By the Envelope Theorem, the ex-ante LMP TTt+i 
at time / -f 1 is given by 

TTt+i = lAt+1 - (2) 

where At+i, gt+i are from ([T]i at time t. 

From we note that the LMP is determined by the 
marginal generator through A and the congestion pattern 
through and /i“. By a congestion pattern we mean the set 
of transmission lines where power flows have reached their 
limits. The forecasting of congestion pattern is a sub-problem 
of the forecasting of LMP. 

III. Multiparametric Programming 

A multiparametric program (MPP) is a mathematical pro¬ 
gram indexed by a vector of parameters. The multiparametric 
programming problem is to solve the mathematical program 
for all values of the parameter vector (or a parameter space 

^By ex-ante LMP we mean that the computation of real-time dispatch and 
associated prices at time t + 1 is based on the estimated system operating 
point at time t and 1-step ahead load and generation forecasts. 


of interest). Consequently, an MPP captures simultaneous 
variations of distributed parameters in a network, which mo¬ 
tivates its application in the problem of LMP forecast. Here 
we summarize some of the key theoretical results essential 
to the development of the proposed probabilistic forecast 
method. The definitions and results follow ca. In particular, 
we are interested in the linear and quadratic programs for 
which the multiparametric programming problems are called 
multiparametric linear programs (MPLP) and multiparametric 
quadratic programs (MPQP) respectively. 

Consider a general right-hand sid^ MPP as follows: 

minz(a:) subject to Ax < b + E6 {y) ( 3 ) 

where x is the decision vector, 9 the parameter vector, z(-) the 
cost functiorfl y the Lagrangian multiplier vector, and A, E, 
b are coefficient matrix/vector with compatible dimensions. 

An MPP solver finds the subset 0 of a given polyhedron 
set such that the mathematical program Q is feasible and 
determines the expression x* (0) G X* (9) of the optimizer (or 
one of the optimizers if it is not unique), where x*(9) is the 
optimal solution and 3C*(9) the set of optimal solutions of Q 
given 9. 

The multiparametric programming analysis builds on the 
concepts of optimal partition and critical region. 

Definition 1 llfTsl'). Let 3 denotes the set of constraint indices 
in (E). Lor any subset J C 3, let Aj and Ej be the corre¬ 
sponding submatrices of A and E, respectively, consisting of 
rows indexed by J. An optimal partition of the index set 3 
associated with parameter 0 G 0 is the partition {3{9),3‘^{9)) 
where 

3{9) ^ G d\A,x*i9) = b + E,9, for all x*{9) G X*( 6 i)}, 
r{9) ^ G 3\AiX*{9) <b + Ei9, for some x*{9) G X*( 6 l)}. 

Lor a given 9o G 0, let (Jo, ^o) — (J(^o), The critical 

region related to the index set Jq is defined as 

03 „ ^ {0 G 0|J(0) = Jo}, 

which is the set of all parameters 0 G 0 with the same active 
constraint set Jq at the optimum(s) of problem (|3ll. 

By Definition [T] the optimal partition specifies two sets of 
constraints by which the congestion pattern is determined: one 
set is a combination of active constraints at the optimum(s) 
and the other inactive. A critical region is the set of all 
parameter vectors for a certain optimal partition. The structure 
of critical region is a key property for the development of the 
online simulation technique, so we present it in Section |V] 
Here, we summarize the global property of critical region and 
its connection with the solution structure of MPP (l3]l. The 
connection between the feasible parameter space and critical 
regions is summarized in the lemma below. 

Lemma 1 (ifTSll. The feasible parameter space 0 can be 
partitioned into critical regions {©i}. If there is no (primal 
or dual) degeneracy, such a partition is unique. 

^By right-hand side we mean the parameter vector 6 is on the right-hand 
side of the constraint inequalities. 

^In this paper, we only consider the cost functions that are linear, piece-wise 
affine or quadratic. 
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Assume that there is no (primal or dual) degeneracy. In the 
linear (or piece-wise affine) cost case, the optimizer is unique 
and the associated Lagrangian multiplier vectors are constant, 
for all parameters within each critical region. Therefore, each 
critical region is associated with a unique LMP vector. In the 
quadratic cost case, the LMP vector is a uniquely defined affine 
function of the parameter within each critical region, which 
follows the theorem below. 

Theorem 1 (ifThl). If the MPQP problem (|5]) is not degen¬ 
erate the Lagrangian multiplier vector associated with the 
optimal solution x*{9) is a uniquely defined affine function of 
6 within each critical region. 

Note that degeneracy may happen for sufficiently large 
systems or high dimensional parameters. Although the La¬ 
grangian multiplier vector may not be unique in such cases, 
a consistent tie breaking rule can be introduced to guarantee 
the uniqueness. 

A key property of the optimizer for both MPLP and MPQP 
is as follows. 

Theorem 2 ( 1151 ifT^ ). If the optimizer x* (9) of MPLP/MPQP 
0 is unique, then it is continuous and piece-wise affine over 
0. In particular, x*{9) is an affine function of 9 in each critical 
region 0^. If there exists multiple optimizers, it is always 
possible to define a continuous and piece-wise affine function 
x*{9) € X*{9) for all 9 G Q. 

This theorem is crucial to the development of the proposed 
forecasting technique. In essence, we treat the randomness as 
a vector of parameter and solve the MPPs ahead of time. The 
online computation of the optimal solution is highly simplified 
because the piece-wise affine mapping between the optimal 
solution and the parameter vector has already been obtained 
offline. 

The problem of solving MPPs, including the computation of 
critical regions and corresponding piece-wise affine functions, 
has been well studied. The first method for solving right- 
hand side MPLPs was proposed by Gal and Nedoma ini. A 
geometric approach for critical region partition was described 
in Ea. The piece-wise affine relationship of the parameter 
vector and the optimal solution was proved for MPLP and 
MPQP in ini and ca, respectively. In this paper, all cal¬ 
culations related to MPPs are performed using MPT3 toolbox 
ifTSll except for dynamical generations of critical regions in the 
proposed online simulation technique presented in Section |V] 

IV. Probabilistic Forecasting Algorithm 

Probabilistic forecasting, in contrast to point forecasting, 
aims to provide the probability distribution of a future LMP. In 
particular, given the estimated system operating point at time 
t and load and generation forecasts, the probabilistic forecast 
at time t of the LMP at time f -f T is given by the conditional 
probability distribution ft+T\t of the LMP vector. Entries of 
the LMP vector are LMPs at individual buses in the system. 
Since each congestion pattern can be mapped from an LMP 

^Note that the cost function of the MPQP problem is assumed to be strictly 
convex. 



Fig. 1. Geometric intuition of the proposed algorithm. 


vector, we only discuss the probabilistic forecasting of LMP 
here; the probability distribution of congestion can be obtained 
from that of LMP. 

The key to probabilistic LMP forecasting is to capture 
spatial and temporal dependencies. Spatial correlations among 
LMPs arise naturally from the optimization that governs the 
real-time dispatch. Temporal correlations, on the other hand, 
are the results of time dependencies in load/generation fore¬ 
casts. The system randomness may also include occurrences 
of random contingencies. In this section, we first give an 
overview of the proposed forecasting algorithm. Details on 
addressing these dependencies are then discussed. 


A. Overview 

The basic idea of the proposed probabilistic forecasting 
technique is using multiparametric programming analysis to 
characterize the variation of the real-time LMP with respect 
to the random load and generation. By formulating the DC- 
OPF problem ([T]i as an MPP in the form of 0, the real-time 
LMP can be expressed as a function of the random load and 
generation. The distribution of the LMP vector at a future 
time can be obtained from the probabilistic forecasts of the 
stochastic load and generation. 

As illustrated in Fig. [T] load and stochastic generation fore¬ 
casts in O are treated as parameters, denoted by 9 = {d, —w), 
where d is the vector of parametric load and w the vector of 
parametric generation (which is treated as a negative load). 
The feasible parameter space 0 is partitioned into critical 
regions {0i,--- , 07 }. Within each region 0^, the optimal 
dispatch is associated with the same Lagrange multipliers, and 
hence a unique LMP vector for all 9 G Qi- Given the 
network parameters, the MPLP solver computes the partition 
{0i}. Correspondences {0i,7ri} are then obtained by the 
Lagrange multipliers and the LMP model 0. Note that this 
computation does not depend on the actual realization of load 
and generation. Therefore, the computation of the partition and 
the correspondences can be obtained offline. 

Consider now the trajectory of a realization of the random 
load and generation process 9t as illustrated in Fig. [T] Given 
the realization 9t and system measurements at time t, we are 
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interested in the conditional probability distributior@ 

ft+T\t{i) = ]P[^i+T s (4) 

As depicted by the shaded circles in Fig. [T] uncertainties 
associated with load and generation forecasts increase with 
time. At time t, the realization of parameter 9t is known 
and thus the distribution ft is an unit vector. But parameter 
9t^T may take values from several critical regions where LMP 
values and congestion patterns are different. Therefore, the 
forecast probability mass function ft+T\t of St+T\t may have 
several non-zero elements. 

The proposed probabilistic forecasting algorithm involves 
two parts explained in the following subsections: the com¬ 
putation of critical regions and the estimation of conditional 
probability distributions, where the former is computed offline 
and the latter online. 


B. Forecast with Varying Operational Conditions 

We describe here a baseline formulation from which critical 
regions are obtained. We formulate the DC-OPF ([T]i used to 
compute LMP at time t + T as the following right-hand side 
MPLP with the uncertainty parameter 9 consisting of only 
stochastic load and generation. 


subject to 


IT 

-IT 


/n 

-In, 







r 0 1 



r 
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IT 1 
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0 

T + 
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1 
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— L~ 

-1 

— St+T-1 

0 

— St + T-1 

0 

^t+T-l 

^t + T-l 



0 

0 


. ~^t + T-l . 


(5) 


where S'i+T-i, and are 

shift factor matrix, vectors of max and min transmission limits, 
and vectors of max and min generation capacities at time t + 
T — 1, respectively. 

Here we allow time varying but known system parameters 
such as shift factors, flow limits on transmission lines, and 
limits on generations, restricting uncertainties only to load 
and generation. The idea is to include deterministically sched¬ 
uled events in the forecasting problem. Examples include the 
scheduled changes in network topology ESI, and generation 
capacity and transmission limit 1201. 


C. Forecast in the Presence of Probabilistic Contingencies 

The baseline MPLP formulation described in Section IIV-BI 
can be extended to include the presence of probabilistic con¬ 
tingencies for unexpected events. For example, a transmission 
line may be tripped in a storm or generation capacity reduced 
due to faults. Such uncertainties in the system parameters 
need to be handled differently from those associated with the 
stochastic load and generation. 

®This is the case for linear or piece-wise affine cost functions. If the cost 
function in the real-time economic dispatch is quadratic, the forecast distribu¬ 
tion ft+T\t of LMP TTt+T can be obtained from the (continuous) distribution 
of 9 and the uniquely defined affine function of the Lagrangian multiplier 
vector within each critical region. 


Because unexpected changes of system configurations are 
typically small probability events, we assume that there are a 
total of K possible system configurations at time t + T. From 
historical data, we assume that system configuration k happens 
with an estimated probability pk- 

We solve the baseline MPLP (|5ll for each system config¬ 
uration and obtain critical regions for system configuration 
k denoted by {0^^^}. By the total probability theorem, the 
probabilistic forecast of LMP at time t + T is therefore given 

by K 

ft+T\t = J2P>^ft+T\V ( 6 ) 

k=l 

(k) 

where *^be forecast distribution under system config¬ 

uration k using critical regions {0j-^^} for fc = 1, 2, • • • ,K. 

We illustrate the above idea with an example. Consider the 
case when at most one of K contingencies can occur between 
time t and t + T. We model the probabilistic contingency 
as independent tosses of a AT -f 1 faced dice where con¬ 
tingency k occurs with probability pk and no contingency 
with po = 1 — further assume that, once a 

particular contingency occurs, it remains until time t + T. We 
then solve each MPLP problem (|5]l for all AT -f 1 possible 
system configurations, including the normal condition. The 
probabilistic forecast of LMP at time t + T is given by 

if contingency k occurs 
Ilk=oPkfl+T\t otherwise 

(k) 

where ft^fu is the forecast distribution under system conflg- 

(k) 

uration k with the feasible space for fc = 0,1, • • ■ , A". 

Note that system configuration 0 denotes the normal condition. 

Having obtained the critical regions, we now consider the 
problem of computing the conditional distribution of LMP at 
time t + T, given the load and generation forecast at time t. 

D. Probability Distribution Estimation 

The estimation of conditional probabilities in (|4]i depends 
on statistical models of load and generation. Such models can 
be obtained from either models of the load and stochastic 
generation process or specific prediction methods used to 
generate load and stochastic generation forecasts. 

As illustrations, we present a directional Gaussian random 
walk model and an autoregressive (AR) noise model for the 
random load and generation processes here. It should be noted 
that any statistical model or prediction method can be applied. 
The purpose of using these models is to gain insights into the 
behavior of forecasting performance by taking advantage of 
some of analytically tractable properties. 

1) A Directional Random Walk Model: We first consider 
a random walk model of the load/stochastic generation based 
on a given mean trajectory. Such a model represents a case of 
minimally informative forecast. Note that the mean trajectory 
can be any available forecast. For example, a reasonable mean 
trajectory is the day-ahead load forecast. 

Assume that load/stochastic generation 9t follows a random 
walk process with a (known) mean trajectory dp. 

= St-i + 9t — 9t-i + et, 



( 8 ) 
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where ct ^ A/"(0, E). Given the realization 9t at time t, the 
actual load/generation at time i + T is given by 

t+T 

dt+T = 9t+ Ot+T — dt + 'Y2 
i— 1+1 

Therefore, the distribution of Ot-\-T conditioning on 9t is; 

AA(0t+T|t,ST) (10) 

where = 9t + 9t+T — St is the conditional mean of 9t+T 

and St = TY, is the cumulative variance within prediction 
horizon T. 

2) An AR Noise Model: The second model we consider 
is an AR(1) noise model where we assume the deviation of 
the load or generation from the expected value is an AR(1) 
process. This is a case when the load or stochastic generation 
is highly structured. In particular, 

9t = 9t + at, at = + et, ( 11 ) 

where 9t is the (known) mean trajectory, (p the parameter of 
the AR process, and et ~ Af(0, E). Given the realization 9t = 
9t + at at time t, the noise at time f + T is given by 

T-l 

at+T = {dt — dt) + (p^'et+T-i, ( 12 ) 

i=0 

and the actual load/generation at time f + T is given by 

T-l 

9t+T = dt+T + {dt — 8t) + 'Y2 ( 13 ) 

i=0 

Therefore, the distribution of 9t+T conditioning on 9t is; 


9t+T ^ J^{&t+T\t,^T) (14) 


where 0 i+T|t = Ot+T + {&t —Ot) is the conditional mean of 

9t+T, and Et = the cumulative variance within 

prediction horizon T. 

To sum up, for both models, the conditional probability of 
9t+T falling in critical region 0 ^ given 9t is: 

exp{-2(x - 9t+T\tyY:p^{x - 9t+T\t)} 
lOi y^(27r)"|ET| 

(15) 


/t+T|t(0 — / 

J0i 


where 0 t+T|t and Et are model associated statistics given in 
(doll and dllli. 

In general, Monte Carlo techniques are necessary to estimate 
the conditional probability (115b . To accelerate the sampling 
process, importance sampling technique is used. In particular, 
for each critical region 0 ^, instead of drawing values from dis¬ 
tribution A/’( 0 t+T|t, Et), we use Af{v{Qi),YT) where u( 0 i) 
is the mean of all vertices of critical region 0^ The estimate 
of /t+T|t(0 is then given by 


N 


ft+T\t(i) ^ ^ y 


i=i 


h{Sj) 


(16) 


where samples {si, • • • , sjv} are drawn from J\f{v{Qi), Et), 
and g{-) and h{-) are probability density functions (PDFs) of 
distribution JV{9^^J•^^,YT) and A/"(u(0i), Et), respectively. 
Note that the importance distribution A/"(v(0i), Et) only 
shifts the mean of the nominal distribution J\f {9 t+T\t, ^ t ) but 
keeps the variance the same. 


3) A Quadratic Cost Case: If the cost function in the DC- 
OPF ([T]i is quadratic, the distribution of the LMP tt^+t cannot 
be estimated by the conditional probabilities of 9t+T falling 
in each critical region. Because the LMP in this case is an 
affine function of the parameter vector by Theorem [T] 

Here we derive the conditional distribution ft+T\t at time t 
of the future LMP itt+T at time t + T given the conditional 
distribution of 9t+T- By Theorem [T] and the LMP formulation 
(| 2 ]i, for each critical region 0 ^, there exists an affine function 
ttii') ■ 0 i —> Hi such that 

TTi{9) = Ui9 + Vi, for all 9 G 0i, 

where H^ is the codomain and Ui and vt are associated 
coefficient matrix/vector. Note that these coefficients can be 
obtained from the affine function of the Lagrangian multiplier 
vector and the LMP formulation (|2]i. 

Given the conditional distribution J\f{9t+T\t, ^t) of 9t+T, 
the conditional PDF of tt^+t is given by 

N 

ft+T\t{Tt) = X]^ni(7r)/t+T|t.ni(7!‘), (17) 

i=l 

where N is the number of critical regions, and ft+T\t.ni ( ) the 
PDF of M{Ui9t+T\t + Vi, UjYTUi), which is the conditional 
distribution of itt+T in codomain H^. 

In summary, the proposed algorithm treats load/generation 
as parameters, formulates the DC-OPF ([T]i as an MPP (|3ll, 
determines the critical regions, and computes the conditional 
distribution of LMP and congestion using load/generation 
forecasts and the real-time operation conditions. These system 
conditions, such as transmission rate, generation capacity, and 
network topology, are allowed to vary with time but known to 
the operator at the time of forecast. Contingency models are 
also incorporated in the proposed technique. 

V. Online Forecasting via Dynamic Critical 
Region Generation 

A limiting factor in the proposed technique above is the 
computational cost associated with the multiparametric pro¬ 
gramming. Since the solution structure is characterized by 
critical regions, all critical regions that partition the parameter 
space have to be calculated. Although such computation can 
be made offline, it may not be computationally tractable for 
large systems because the number of critical regions may grow 
exponentially with the number of constraints. 

In this section, we propose an online Monte Carlo technique, 
referred to as dynamic critical region generation (DCRG), 
that significantly reduces the computation cost. The idea is 
to take advantage of the fact that, in practice, random load 
and generation processes are bounded and tend to concentrate 
around their mean trajectories. As a result, a small fraction 
of critical regions represents the overwhelming majority of 
observed critical regions. When a parameter falls in a critical 
region that was visited before, the Lagrange multipliers that 
are used to generate LMPs can be obtained directly from the 
affine mappings associated with that critical region, without 
solving the DC-OPF ([TJ. 
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The key idea of DCRG, therefore, is to compute on demand 
the critical region and the associated coefficients of the affine 
mapping of parameter to the LMP. This computation, fortu¬ 
nately, is no more than elementary matrix inversions and mul¬ 
tiplications. The computation procedure is given by the follow¬ 
ing theorem that summarizes known results in MPLP/MPQP 
na Qa. The proof is provided in the Appendix. 

Theorem 3. Given a parameter Oq, let {3 q,3q) be the optimal 
partition of the index set 3 of (0 associated with 9 q. Let 
Ajg , Ejg and bjg be, respectively, the submatrices of A, E 
and subvector of b corresponding to the index set Jq- 
Ajc , Ejc and bjc be similarly defined for the index set Jq. 
Assume that MPP 0 is neither primal nor dual degenerate 
for all 0. Denote the critical region that contains 9 q by © 3 ^,. 

(1) For the MPLP 0 with cost function z(x) = c'^x, the 
critical region ©j^ is given by 

03o = . (18) 

(2) For the MPQP 0 with cost function z{x) = ^x'^Hx 
where FI is positive definite, the critical region ©j^, is 
given by 

©3„ = (19) 

where Tp and Td are polyhedra defined by 

T A,.H-^A\{A,gH-^A\)-\hg+E,g0) 

Td = {0\iAjgH-^Al)-^\bjg +"Ejg0) ^ O}. 

In applying DCRG to LMP forecasting using online Monte 
Carlo simulation, we generate samples of random load or 
generation, either based on load/generation models or from 
historical data. Instead of directly simulating the real-time 
market operation as in 0 , we check if the generated sample 
falls into a critical region that has been used before. If it 
does, the LMP can be generated using the affine mapping of 
Lagrangian multipliers directly without solving the DC-OPF 
O. If the parameter does not belong to any critical region 
in the database, a DC-OPF is solved and a critical region 
containing the parameter is computed according to (fTsT l for 
linear cost case or (fT^ for quadratic cost case in Theorem |3] 

VI. Evaluation 

In this section, we present simulation results to compare per¬ 
formances of the proposed probabilistic forecasting algorithm 
with some benchmark methods. We first show results of a 3 
bus system with a linear cost function to gain insights into the 
behavior of the proposed algorithm under various scenarios. 
Simulations using the IEEE 118 bus system with a quadratic 
cost function are then presented to demonstrate the scalability 
of the proposed algorithm and the effectiveness of the online 
heuristic approach given in Section |V] 

A. Benchmarks and Performance Measure 

We compared the proposed techniques with some existing 
benchmarks for forecasting and computation performance. 
Since, to our best knowledge, there is no probabilistic fore¬ 
casting techniques for the operator in the literature, we used 
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the direct Monte Carlo simulation method proposed in 0 
as the probabilistic forecasting benchmarlfl. We also included 
comparisons of the proposed technique with two well known 
point forecasting methods to illustrate the performance gain. In 
particular, the deterministic prediction uses the mean trajectory 
of load/stochastic generation Ot+T to calculate LMP and 
congestion pattern at time t + T. The certainty equivalence 
prediction incorporates measurements at time t and uses the 
conditional mean trajectory 9t+T\t- 

Before presenting numerical results, we introduce a per¬ 
formance evaluation metric of probabilistic forecasts. The 
LMlj^congestion pattern is a discrete random vector. The 
probabilistic forecast of such a random quantity belongs to the 
so-called categorical forecast, and its performance is measured 
by the consistency as well as the statistical concentration of 
the forecast. A standard metric IHI for this type of forecast is 
the Brier Score (BS) 1221 that measures the average distance 
( 2 -norm) between the forecast distribution ftJrT\t and the 
point mass distribution at the realized random variable tti+t- 
Specifically, 


BS(/t+T|t) = E\\ft+T\t - (20) 

where the expectation is taken over all randomness between 
time t and t + T. In (I201 i. ft+T\t is the conditional probability 
vector whose ith entry is given by ft+ritil) = Pi'('Tt+T = T^i), 
and S(x) is the unit vector that is one at entry x and zero 
elsewhere. This score ranges from 0 for a perfect forecast to 
2 for the worst possible forecast. 

Since BS is a succinct formula to measure the overall per¬ 
formance in terms of uncertainty, reliability and resolution, we 
also provide a more intuitive assessment — reliability diagram. 
Reliability diagram is a graph of the observed frequency of an 
event plotted against the forecast probability of an event. It 
measures how closely the forecast probabilities of an event 
correspond to the actual chance of observing the event. 


B. Case Study: A 3 Bus System 

Consider a 3-bus system as depicted in Eig. |2l Generator 
incremental costs and capacity limits are presented in the 
figure. All lines are identical with the maximum capacity of 
100 MW. 

^We want to point out that the direct Monte Carlo simulation approach 
generates exactly the same probabilistic forecast as the proposed technique. 

*Note that the discreteness of LMP is only for linear or piece-wise affine 
cost functions. 




TABLE I 

Critical regions, LMPs and congestion patterns for the 

BASELIN^. 



Critical Region 

LMP 

Congestion 

1 

(0, 130) 

(10, 10, 10) 

(0, 0, 0) 

2 

(130, 170) 

(15, 15, 15) 

(0, 0, 0) 

3 

(170, 200) 

(10, 20, 15) 

(1, 0, 0) 



Time lndex(5-min Interval) 
(a) Random walk model 



Time lndex(5-min Interval) 
(b) AR(1) noise model 


Fig. 3. Impact of load statistical models. 

1) Baseline: We first evaluated the baseline algorithm with 
the two load models described in Section IIV-DI Note that 
only the load at bus 2 was stochastic, which gave the one 
dimensional parametric linear program. 

Table 1 shows the critical regions of the parametric linear 
program, and the associated LMPs and congestion patterns. In 
this case, the parameter load d at bus 2 was partitioned into 
three segments D = {(0,130), (130,170), (170, 200)}. 

We used a straight line ranging from 110 MW to 190 
MW with 2 MW increments as the mean trajectory of load 
dt- The coefficient (f) in AR(1) noise model was set at 0.9. 
The independent noise sequence e* in both models followed 
the standard normal distribution, i.e., A7(0,1). Monte Carlo 
simulations were used to obtain estimated BSs in Fig. [3 

From Fig. [3 we observed that the proposed probabilistic 
forecasting algorithm “Alg-P” consistently outperforms the 
deterministic “Alg-D” and certainty equivalence “Alg-C” pre¬ 
dictors in both load models. The superior performance of the 
proposed technique in these two extreme models (a minimally 
informative model and a highly structured model) shows its 
capability of incorporating different load forecasting methods 
and its forecasting power. Two interesting phenomena are 
worth closer examinations. First, for both load models, peaks 
occurred at the boundaries of two neighboring critical regions 
when the mean load dt is 130 at f = 10 and 170 at f = 30. At 
the boundary point, the probability of dt+r falling in either 
the left or the right critical region was the same. Roughly half 

^In Table |I] the triple for LMP contains price values at bus 1-3. For 
congestion, the triple summarizes status of line 1-2, line 1-3, and line 2- 
3, respectively, where “1” indicates positive congestion (the flow reaches the 
line limit in the positive direction), “—1” negative congestion, and “0” no 
congestion. 



Forecast probability 


Fig. 4. Reliability diagrams for critical region 1 with different load forecast 
uncertainties. 

of the time Alg-D predicted the LMP correctly, the other half 
was completely wrong. From the definition in (120b . the BS of 
Alg-D should be 1. 

Second, a more subtle point, the scores for the random walk 
model showed a slight asymmetry with respect to boundaries 
of neighboring critical regions; the BSs of Alg-P and Alg-C 
at the second peak were lower than that of the first peak, and 
the ranges of non-zero score neighborhood of the second peak 
(from time 24 to 39) of all three algorithms were wider than 
that of the first peak (from time 5 to 16). This phenomenon 
arose primarily from the process of generating the sample 
trajectory of load. Since the entire sample trajectory was 
generated at once, the deviation \dt — dt \ from the mean grows 
over time which leads to a bigger variance of the time crossing 
the second boundary 170 than that of the time crossing the first 
boundary 130. In other words, comparing to the probability of 
crossing the first boundary 130 at f = 10, the probability of 
crossing the second boundary 170 at f = 30 is lower, but the 
probability in its neighborhood is higher. Therefore, the scores 
of Alg-P and Alg-C at the second peak were lower, and the 
ranges of non-zero score neighborhood bigger. 

To evaluate the robustness and reliability of the proposed al¬ 
gorithm, we tested different levels of load forecast uncertainty. 
Specifically, we varied the variance of the noise in the random 
walk load model: by default, e A/"(0, CTq), where ctq = 1, and 
then we used the distribution A7(0, cr^) with different values 
of (T^. The reliability diagrams for the probabilistic forecast 
of critical region 1 using different load forecast uncertainties 
are given in Fig. |4] The results show good resolution at the 
expense of reliability. 

2) Forecast with probabilistic generation outage: We now 
considered the case of a generating unit outage with a partial 
loss of capacity, assuming that the maximum capacity of the 
generator at bus 1 can be reduced to 100 MW with probability 
p. Other settings were the same as those in the first scenario. 

The critical regions under this configuration became D™* = 
{(0,100), (100, 200)}. To predict a future price, we considered 
all critical regions {D, D™*} that load dt+r may fall in, 
where D refers to the critical regions in Table |T] under normal 
conditions. For the outage frequency p, we chose two levels; 
p = 0.01 and p = 0.1. The random walk model was adopted 
to generate load profiles. As benchmarks, both deterministic 
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Time lndex{5-min Interval) 



Fig. 6. LMPs at 10% and 20% renewable penetration levels. 


Fig. 5. Impact of outage frequency p. 


and certainty equivalence forecasting algorithms also took 
contingencies into consideration. 

From Fig. |5] the behaviors of all three algorithms with the 
same outage frequency p were similar to that in the first sce¬ 
nario. Boundary effects and peak asymmetries were observed 
with contingencies as well. Compare the performances of each 
algorithm with different outage frequencies: the smaller the 
outage rate the better the forecast. 

C. Case Study: IEEE 118 Bus System 

The IEEE 118 bus system was used to show the scalability 
of the proposed technique and the effectiveness of the heuristic 
algorithm described in Section |V] The simulations results were 
focused on the complexity comparison between the proposed 
techniques and the probabilistic forecast benchmark. 

We introduced 12 wind generators (roughly 10% of buses) 
with 10 ~ 20% renewable penetratioio The wind generators 
were located at bus 25, 26, 90, 91, 100, 103, 104, 105, 107, 
110, 111, and 112. The selection of these locations were 
intended to represent two wind farms, the small one has 2 
wind generators, and the large one has 10 wind generators 
concentrated on a few neighboring buses. All wind generators 
were assumed to be identical with the maximum capacity of 
110 MW. Denote the wind generation space by the hypercube 
W G We imposed the maximum capacity of 100 MW on 
transmission line 8, 126, and 155. The load profile, generator 
capacities and cost functions, and line and bus labels were 
referred to as in “easel 18” in MATPOWER 1231 . Note that 
the cost function in this system was quadratic, thus the LMP 
was an affine function of the wind production within each 
critical region. 

The mean trajectory w{i) of each wind generator i was 
assumed to be linear, for i = I, - - ,12. In particular, the 
trajectory starts from 10% penetration level, i.e., WQ{i) = 
35.35, at time 0, and ends at 20% penetration level, i.e., 
wioi'i) = 70.70, at time 10. The increment was assumed to 
be constant, i.e., 3.535 MW. LMP values at 10% and 20% 
penetration levels are given in Pig.|6] We observed that higher 
renewable penetration reduced LMP at most buses, but raised 
the LMPs at bus 93, 94, 95, and 96. The reason of such 
nonintuitive increase was the congestion of transmission line 

*®By x% renewable penetration we mean the mean value of total wind 
generation is x% of the total electricity load (4242 MW in this system). 
Note that load was assumed to be deterministic in this case. 
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Fig. 7. Predicted LMP distributions and Gaussian approximations at bus 49, 
90, 94, and 100. 


155 caused by the increased wind production from the big 
wind farm. The random walk model was used for the stochastic 
generation profile. The distribution of the independent noise 
process e* was set to be the standard multivariate Gaussian. 

Results for the 10-step ahead prediction at time t = 0 are 
provided in Pig. [T] The predicted marginal distribution of LMP 
at bus 49 exhibits the a Gaussian distribution as it is well 
fitted by the Gaussian approximation with the sample mean 
and sample variance as distribution parameters. However, such 
Gaussian characteristics were not observed on the other three 
buses. According to the distributions of LMP at bus 94 and 
100, the extreme values of LMP occasionally appeared as 
spikes, which were caused by the network congestions. 

In the following, we will focus on the complexity perfor¬ 
mance of the proposed forecasting techniques. Theoretically, 
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Critical Region Index 

Fig. 8. Conditional distribution of observed critical regions at time 10. 


TABLE II 

Computation TIME (in seconds) for 10,000 samples. 



Offline 

Online 

Total 

Alg-P 

160.60 

23.32 

183.92 

Alg-DCRG 

— 

23.59 

23.59 

Alg-MC 

— 

52022.57 

52022.57 


computation, demonstrating the efficiency for the real-time 
LMP forecasting. 



Fig. 9. The expected number of DC-OPF computations vs. the total number 
of Monte Carlo simulations. 


there were at most 2^^® critical regions, because the constraints 
in the DC-OPF ([TJ for this example included 1 energy balance 
constraint, 6 transmission constraints, and 108 generator ca¬ 
pacity constraints. For the given hypercube W of the parameter 
space, there were in total 273 critical regions. But for the 
generated 10,000 samples at time 10, only 17 critical regions 
were observed and more than 99% samples fell in 3 critical 
regions, as shown in the distribution over the observed 17 
critical regions in Fig. |8] 

Instead of exploring the entire wind production space of¬ 
fline, we implemented the online forecasting algorithm DCRG 
given in Section IV] Fig. |9] shows the comparison of computa¬ 
tional cost between the proposed DCRG algorithm, labeled 
by “Alg-DCRG”, and the direct Monte Carlo simulations, 
labeled by “ Alg-MC”. Both algorithms present approximately 
linear growth in the logarithm scales. But the proposed DCRG 
algorithm provided more than three orders of magnitude 
reduction in the number of DC-OPF computations required 
in the simulation. 

Finally, we provide the computation time comparison for 
the three probabilistic forecasting techniques in Table HH 
All computational times were evaluated by implementing the 
algorithms in Matlab environment with the default “quadprog” 
solver and an external MPT3 ca toolbox on a desktop 
with an Intel Core 17-3770 CPU at 3.4 GHz and 8 GB 
memory. No attempts were made to optimize the efficiency 
of the algorithms and their simulations. From Table HH we 
can conclude that the direct Monte Carlo simulation approach 
Alg-MC failed to meet the time constraint for real-time LMP 
forecasting, since 10,000 samples took more than 14 hours 
to generate the distribution. The proposed techniques, on the 
other hand, only took less than half a minute in the online 


VII. Conclusion 

This paper presents a new methodology for the short-term 
forecasting of real-time LMP and congestion for system op¬ 
erators. Based on a multiparametric programming formulation 
of DC-OPF, we have developed an approach that exploits the 
parametric structure of DC-OPF solutions to obtain conditional 
distributions of future LMP and congestion. 

For system operators, the proposed online forecasting tech¬ 
nique provides a new tool for managing operation risks and 
solving stochastic optimization problems. See, for example, 
0 . For market participants, on the other hand, congestion and 
LMP forecasts by the operator provide actionable signal for 
managing flexible resources and demand side management. 

Currently, there are very few probabilistic forecasting tech¬ 
niques for system operators beside direct Monte Carlo simu¬ 
lations JS). The approach presented here represents a first step 
toward online forecasting in large power systems with signif¬ 
icant stochastic components. For future work, there is a need 
to develop computationally tractable techniques, informative 
performance measure, and a set of practical benchmarks. 

Appendix 

Proof of Theorem[3] 

Proof. For any parameter 9, denote the optimal primal and 
dual solutions to MPP 0 by x*(9) and y*(9) respectively. 

(1) For the MPLP case, if 9 is in the same critical region as 
9o, they have the same optimal partitions, which means that 

Aj^x*{9)-bj^-Ej^9 = 0 , ( 21 ) 

Ajcx*{9)-bjc-Ejc9 < 0 . ( 22 ) 

Because MPLP is neither primal nor dual degenerate, Aj^ has 
full rank, and 

x*{9)=Af^\bj„+Ej,9). (23) 

Substituting x*{9) into (122b . we have 9 G ©Jq. 

(2) For the MPQP case, the first order optimality conditions 
are given by 


Hx*{9)+A^y*{9) = 0, (24) 

Ax*{9) - b- E9 < 0, (25) 

y*{9){A,x*{9)-h-E,9) = 0, Vi € (26) 

y*{9) > 0. (27) 

From (124b . 

x*{9) = -H-^A^y*{9). (28) 
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Substituting the result into (|2^ , we have 

A,H-^AjyU0)+h + E,e = 0, VzGJo, (29) 

y*ie) = 0, Vi G Jg. (30) 

By the non-degeneracy assumption, the rows of Aj^ are 
linearly independent. This implies that AjgH~^Aj^ is a square 
full rank matrix. Therefore, from ( I29l l, we solve 

y* 3 „{S) = + E,,e) (31) 

and substitute y^ (0) and yjc (0) into (l28l l to obtain 

x*ie) = H-^Al{Aj,H-^Al)-\bj,+Ej,9). (32) 

Substituting x*{6) from (132b in the primal feasibility con¬ 
ditions (l25l l gives Tp and substituting x*{9) from (l32T i in the 
dual feasibility condition (|27] | gives Td- We therefore have 
9 G 0jo. □ 
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